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Abstract 

Ground state properties and phonon dispersion curves of a classical linear 
chain model describing a crystal with an incommensurate phase are stud- 
ied. This model is the DIFFOUR (discrete frustrated (j) 4 ) model with an 
extra fourth-order term added to it. The incommensurability in these models 
may arise if there is frustration between nearest-neighbor and next-nearest- 
neighbor interactions. We discuss the effect of the additional term on the 
phonon branches and phase diagram of the DIFFOUR model. We find some 
features not present in the DIFFOUR model such as the renormalization of 
the nearest-neighbor coupling. Furthermore the ratio between the slopes of 
the soft phonon mode in the ferroelectric and paraelectric phase can take on 
values different from —2. Temperature dependences of the parameters in the 
model are different above and below the paraelectric transition, in contrast 
with the assumptions made in Landau theory. In the continuum limit this 
model reduces to the Landau free energy expansion for type II incommen- 
surate crystals and it can be seen as the lowest-order generalization of the 
simplest Lifshitz-point model. Part of the numerical calculations have been 
done by an adaption of the Effective Potential Method, originally used for 
models with nearest-neighbor interaction, to models with also next-nearest- 
neighbor interactions. 

PACS numbers: 64.70.Rh, 64.60.Kw, 02.60.Pn 
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I. INTRODUCTION 



Like other transitions, the phase transition to an incommensurate (INC) phase (reviews 



essentially in accounting for the expansion of the free energy density as a function not only 
of the components of the order parameter, but also of their spatial derivatives. Therefore 
the global free energy becomes a functional of spatially dependent components of the order 
parameter and the equilibrium configuration for given values of temperature and external 
parameters is found as a solution of a variational problem. 

The continuum Landau theory allows a natural classification! of the possible forms of the 
free-energy functional for an INC transition into two classes, according to whether the driving 
term in the free energy expansion responsible for the appearance of the incommensurate 
state is linear (type I, Lifshitz invariant present) or quadratic (type II, no Lifshitz invariant) 
in the gradient of the order parameter. The properties of those two kinds of INC phases 
are different: for type I INC phases the lock-in transition is either continuous, or only 
slightly discontinuous, and approaching the lock-in temperature T c these phases exhibit the 
structuration of the modulated phase into discommensurations or solitons. On the other 
hand, the modulation of the type II INC phase remains practically sinusoidal down to the 
temperature T c and the lock-in transition is always of first order. Although the above 
statements can be considered as a rule of thumb, there are cases known where there is 
coexistence of solitonic and sinusoidal structural modulation. See Aramburui for details. 

In the following we will only be concerned with models describing type II INC phases. 
Landau theory has been rather successful in describing basic properties of these phases, but 
if one wants to have a better understanding of the true microscopic origin of the INC phase 
one has to go beyond this phenomenological approach. One possibility would be to study full 
microscopic models with realistic interactions. Another approach, to which the main part 
of this paper will be devoted, is to study semi-microscopic models which take into account 
the discrete nature of the systems and discuss properties in terms of (effective) inter-atomic 
interactions. 

The discreteness of a lattice leads to a number of important physical consequences such 
as pinning of solitons in the anisotropic next-nearest-neighbor Ising (ANNNI) model (see 
YeomansEi] for a review) and in the Frenkel-Kontorova modeljl and the existence of a devil's 
staircase (infinite number of commensurate and incommensurate phases), for example found 
in betaine calcium chloride dihydrate (BCCD)i Within Landau theory it is difficult to 
explain the occurrence of a specific sequence of transitions: several lock-in terms are needed 
then. For example, Ribeiro et a/0 chose the magnitude of four distinct lock-in contributions 
to the free energy in such a way as to stabilize the four most prominent commensurate phases 
in BCCD. Furthermore, in discrete models chaotic states are possible,00 which may provide 
an alternative description of phenomena observed in for example spin glasses, superionic 
conductors, the magnetic system CeSb and systems with pinning of charge density waves. 
See Bak§ for details. 

In the past few years different lattice models have been constructed, for example to 



crystals with Pcmn symmetryllil These models are 2-dimensional, with only nearest-neighbor 





BCCD§ and, more general, in 
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interactions. Hlinka et alx3 studied a 3- dimensional nearest-neighbor model, applicable to 
BCCD. All these models have in common that the frustrated interaction, needed for having 
an incommensurate phase, comes from a nearest-neighbor mixing interaction. 

Recent X-ray,@ neutron0lli and Raman0 experiments on the Sn 2 P2(Si_ x Se x )6 crystal 
family of uniaxial ferroelectrics motivated us to study lattice models. In the composition- 
temperature phase diagram of this crystal family a Lifshitz point is present.il At this point 
the paraelectric phase, the ferroelectric phase and the INC phase become equal, and the 
boundaries separating these phases have equal derivative. Such a special point was, except 
for some ferroelectric liquid crystals, only found in the temperature-applied magnetic field 
phase diagram of the magnetic compound MnP.ll] From both experimental and theoretical 
point of view this uniaxial Lifshitz point is interesting because critical exponents deviate!! 
from those found for ordinary critical points. This crystal family furthermore displays an 
interesting modulation wave vector behavior, shows cross-over effects from order-disorder to 
a displacive type of phase transition,il and the ratio between the slopes of the soft phonon 
mode in the ferroelectric and paraelectric phase deviates^ substantially from the standard 
value R = —2. These features are much easier to understand in a lattice model than in 
Landau theory. 

The paper is arranged as follows: in Section II we present a one- dimensional model; in 
many anisotropic systems, like Sn 2 P2(Si_ I .Se a: ) 6 , the modulation wave vector is in one specific 
direction. The incommensurability may arise if there is frustrating interaction between 
nearest-neighbor and next-nearest-neighbor couplings. We discuss general features of the 
model. In Section III we give some exact results regarding ground state properties. In 
Section IV we discuss the dynamics and the stability of the various phases. In Section 
V the phase diagrams, calculated partly analytically, partly numerically, are presented. 
Temperature effects are treated in Section VI. In Section VII we discuss the continuum limit 
of the model which gives the connection with the Landau theory. We conclude and give 
an outlook for further research in Section VIII. In Appendix A we give some exact results 
for phase boundaries and in Appendix B we present the next-nearest-neighbor extension of 
the Effective Potential Method for the determination of the ground state. This method was 
used to calculate some of the phase diagrams. 



II. AN EXTENSION OF THE DIFFOUR MODEL 

In the following we will be concerned with an extension of the so-called DIFFOUR 
modelfEl (discrete frustrated (f) 4 model), for which the potential energy can be written as 

f 2 B a C . s 2 D . .2 

n 

E 1 

+ Y [ x l ( X n ~ Xn-lf + xl (X n - X n+ if] j. (2.1) 

The original DIFFOUR model, or EHM (elastically hinged molecule) model,i! has E = 0. 
Although in principle this model gives incommensurate ground states, the behavior of the 
modulation wave vector as found in experiments can not be reproduced satisfactorily by the 
model. In order to account for this shortcoming we supply the DIFFOUR model with a 
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non-linear coupling to neighbors. There are several possibilities: if we restrict ourselves to 
fourth-order terms we can consider a term oc (x n — x n _i) 4 . The resulting model has been 
studied by Lamb,@ who showed that the origin of this term is related to strain terms in 
the so-called magnetoelastic DIFFOUR model. Another possibility would be to consider a 
term of the form ) . Without the term mentioned above this would however 

be rather unphysical. Instead we will choose a term of the form oc x 2 n (x n — x n ±i) 2 . This is 
the lowest-order dispersive fourth-order term, as will be shown in Section VII, and can, for 
example, be obtained from strain terms. 

The order parameter x n can be, for example, a displacement, a component of the po- 
larization P for ferroelectric systems, a component of the magnetization M for magnetic 
systems, a rotation angle or a strain component. In this article we use for convenience terms 
like paraelectric, ferroelectric and antiferroelectric to distinguish between different ground 
states. The origin of incommensurability in this model is essentially competition between 
interactions with nearest- and next-nearest neighbors which may lead to frustration. Higher 
order terms are needed for stabilization. 

We expect that the extra fourth-order term (E ^ 0) has a large effect on the phase 
diagrams for E — 0. To actually determine phase diagrams it is not necessary to vary 
all 5 parameters A, B, C, D, E, which can be seen as follows: by taking x' n = y/B/\D\x n 
and V = B/\D\ 2 V we get the following renormalized parameters: A' = A/\D\, B' = 1, 
C' = C/\D\, D' = D/\D\ = ±1, and E' = E/B. 

For some purposes it is convenient to rewrite the potential in the following form: 



r A 2 B 4 

^ — x n -\- —x n -\- C X n X n ^i -\- Dx n X n —2 



E ~i 

+ Y K ( X n ~ x n-lf + xl (x n ~ X n+1 f] j, (2.2) 

with A = A + 2C + 2D, B = B, C = -C, D = -D and E — E. From this the connection 
with the ANNNI model can easily be made. Let us put E = and B = —A. If we now 
take the limit A ^ -oo we end up with a model with two infinitely deep wells. The x n can 
only take on values ±1 and can thus be seen as spins. These spins are coupled to nearest- 
neighbors and next-nearest-neighbors via the C and D terms. So by increasing the depth 
of the double-well potential there is a crossover from displacive behavior to order-disorder 
behavior in the transition from the normal to the incommensurate phase. 
Inserting x' n = (— l) n x n in the above potential leads to 



^ = S ^/\2 X n + ~^ X n ~ ^ X n X n-l + ^ >X n X n-2 



E 

+ 2 



X 



n { X n + X n-l) + X n ( X n + X n+l) \- (2-3) 



In the DIFFOUR model (E = 0) this leads to the following symmetry: if {x n } is a state for 
C = X, then {(— l) n x n } is a state for C = —A with the same energy. This property can 
for example be seen in the ferroelectric-antiferroelectric phases. However, for E 7^ this 
symmetry C <-> — C is no longer present. 
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III. GROUND STATE PROPERTIES 



Different ground states are possible, depending on the values of the parameters. The 
stationary states are solutions of dV/dx n = 0, giving 

Ax n + Bx\ + C {2x n — x n _i — x n +i) + D (2x n — x n _2 — x n+2) 

+E [Ax\ - 3x 2 n (sn-i + ^n+i) + 2x n (x 2 n _ x + - (x^ + = 0. (3.1) 

If we impose periodic boundary conditions xn+u = %n we arrive at a set of N coupled non- 
linear equations. To find the lowest-energy state for each solution and for each value of N 
the potential energy has to be evaluated. For low-period commensurate states (small values 
of N) analytic solutions of the above equation can be found. In the following, we study (for 
fixed N) periodic solutions of ( |3.1|) . For them we give the equilibrium values {x n } and the 
corresponding energy per particle v = V/N. 

In the paraelectric state (N = 1) all particles are in the equilibrium positions 

x n = 0, v = 0. (3.2) 

In the ferroelectric state (N = 1) all particles are uniformly displaced from their equi- 
librium positions 



x »=V"f v = ~^- (3 - 3) 

Note that B always has to be positive, for the potential to be bounded from below. This 
implies that the ferroelectric state only can exist for A < 0. For A > the ground state 
may be paraelectric. 

In the following we give some analytic results, based on numerical calculations of the 
shape of the solution {x n }. 

In the antiferroelectric state (N = 2) particle positions alternate along the chain 



A + iC (A + iCf 

Xn = ' ' ^ ~B~+WE' P= - 4(B + 16E) - (3 ' 4) 

The potential is unbounded from below for E < —5/16 and stable solutions exist only for 
E > —B/16 and A + AC < 0. Both conditions have to be satisfied. 

For iV = 3 we determined the solution with lowest energy to be of the form (xi,x%, £3) = 
with X\jx-i < 0, and 

2 = A + 2(C + D)(l-k) _ Ak + (C + D)(k-1) 
5 " B + 2E{2 - 3k + 2k 2 - k 3 ) Bk 3 + E(2k 3 - 3k 2 + 2k - 1) ' 1 ' ' 

The factor k is determined by 

2 [B(C + D) + E(-A + C + D)]k A - [B(A + 2C + 2D) + 2E(-A + 2C + 2D)] k 3 (3.6) 
—3AEk 2 + [B(A + C + D) + 2E(A + 2C + 2D)} k - [B(C + D) + E(-A + 2C + 2D)] = 0. 

The energy per particle is given in terms of £ and k as 
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v = f U(l + 2k 2 ) + 2(C + D)(l - A;) 2 ] + f- [5(1 + 2A; 4 ) + 4£(1 + k 2 )(l - kf] . (3.7) 
6 12 

Note that the quartic equation (|3.6| ) can be written as 

(1 - Jfe){-2 [B(C + D) + E(-A + C + D)}k 3 + [BA + 2E[C + D)] k 2 

+ [BA + E(3A + 2C + 2D)] k - [B(C + D) + E(-A + 2C + 2D)} } = 0, (3.8) 

where the special solution k — 1 gives a ferroelectric state. The solution of the remaining 
cubic equation, which can be solved exactly for given parameters, gives a k such that X\jx^ < 
0, a true N = 3 state. 

The lowest energy state for N = 4 has X\ = X2 = p, £3 = £4 = — p with 



, A + 2C + 4D (A + 2C + AD) 2 
P= \l B + 8E ' W = 4(5 + 8E) • (3 ' 9) 

We have to keep in mind that we must satisfy E > —5/16, which is not obvious from the 
above expression, but comes from the analysis of the antiferroelectric state. 

The lowest energy solution for N = 6 can analytically be obtained, in the same manner as 
for N = 3. It has the form (xi, X2, X3, £4, x$, Xq) = (k£, £, k£, —k!;, —k£). The lowest en- 
ergy solution for iV = 8 reads (xi,X2, £3, £4, £5, xq, X7, xs) = (k£, k£, —k£, —k£). 
For N = 5 one needs 2 different values of k: (xi, x 2 , x 3 , x 4 , x 5 ) = (k£, k'£, £, k'£, k£) and for 
N = 7 one needs three different values: k, k', k", and the above method no longer works for 
N = 5 and N = 7. Therefore, to find states with larger periods, or even incommensurate 
periods, we rely on numerical calculations, for which true incommensurate states of course 
never can be found. However, the idea is that such a state can always be arbitrary well 
approximated by a commensurate state with wavelength 

N 

A = — , s = l,2,..., iV>2s, (3.10) 
s 

where N, s are coprime numbers. Such a solution has a period N and, in general, 2s = 
(number of local minima + number of local maxima). In the special case where the {x n } 
take on positive and negative values, 2s = (number of sign changes within the period N). 
The bigger N and s, the better the approximation. 

As an example of a numerical calculation we consider the ground state for A = 
2.24999, B — 1, (7 = 1, D = —1, E — 1. It is known to be incommensurate (see Sections 
IV and V) with a wavelength arccos(|) ps 4.76679213. By the Farey construction! we find 
that || ~ 4.76923077 should be a reasonable commensurate approximation. We numerically 
determined the ground state in terms of the {x n }. The result is shown in Figure |l]. After 
62 particles the sequence repeats itself, and the solution passes through zero 26 times. We 
can label this state by its modulation wave vector ||, measured in units of 2-7T. 

For certain regimes in the parameter space the ground state can be determined analyti- 
cally. The entire phase boundary of the paraelectric state and a part of the phase boundary 
of the ferroelectric state can be calculated. For the nearest-neighbor case in the DIFFOUR 
model proofs are given by Janssen and Tjon.0 We extend their proofs to the case in which 
we also have next-nearest-neighbor interaction. As the proof is rather lengthy it will be 
given in Appendix A. 
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IV. PHONON DISPERSION CURVES AND STABILITY LIMITS 



To decide whether a solution of the equilibrium conditions is locally stable or not, one 
considers small displacements e n from the positions given by a static solution {x n } satisfying 

u n = x n + e n . (4.1) 

The phonon frequencies are given by the square roots of the eigenvalues of the dynamical 
matrix and for stability all eigenvalues have to be non-negative. The dynamical matrix for 
a period-A solution (N > 5) has elements 

D njn = A + 3Bx 2 n + 2C + 2D + 2E [6x 2 n - 3x n (x n ^ + x n+1 ) + x 2 n _ x + x 2 n+1 ] , 
D n>n±1 = -C + E [-3x 2 n + Ax n x n±1 - 3x 2 n±1 ] , (4.2) 

D n ,n±2 — —D, 

with xn+u = x n and the special cases 

D ljN — [-C + E {-3x\ + Ax x x N - 3x 2 N )] e~ iq , 
D^n-x = D %N = -De~ iq . (4.3) 

In the above expressions the {x n } are solutions of ( |3.1| ). Furthermore -D n , m = D* mn and all 
other matrix elements are zero. 

In the following the phonon branches for certain low-period states will be examined. This 
will be done in terms of A, B, C, D, E. 



A. Paraelectric state 

For the paraelectric state (|3.2|) the dynamical matrix is given by 

D = A + 2C(1 - cos(g)) + 2D(1 - cos(2g)) = mu 2 . (4.4) 

Rewriting this equation gives 

muj 2 = A + (AC + 1QD) sin 2 (|) - 16D sin 4 (|). (4.5) 

Note that this expression does not contain E. This means that the stability limits for 
the paraelectric phase in this model are the same as those for the paraelectric state in 
the DIFFOUR model. Now, we are looking for the minimum of this phonon branch. We 
distinguish the cases D > and D < 0. The results are summarized in Table |. 

For D > and C = 0, the branch has two minima, at sin 2 (|) = and sin 2 (|) = 1. For 
A < the ferroelectric state and the antiferroelectric state are degenerate, for E = only. 
Comparison with calculations in Appendix A shows that as long as the paraelectric state is 
stable, it is the ground state. Destabilization is the condensation of a soft phonon. 



7 



B. Ferroelectric state 



For the ferroelectric state (13.31) one has 



muo 2 = A + 3Bx 2 + 2C(1 - cos(g)) + 2D(1 - cos(2g)) + AEx 2 (l - cos(g)) 

= -2A + {AC + IQD - — — ) sin 2 (^) - 16L>sin 4 (^). (4.6) 

B 2 2 

See Table [Tl] for the analysis. Note that for D > there is degeneracy for E = 0. 

C. Antiferroelectric state 

Finally, the phonon branches of the antiferroelectric state ( |3.4j ) will be investigated. The 
2x2 dynamical matrix is given by 

A 4- AC 

D ltl = L>2, 2 = A + 2C + 2D(1 - cos(g)) - (35 + 28E)- 



B + 16E 1 
A 4- AC 

Di,2 = D* 2>1 = {e* + l)(-C + WE^^). (4.7) 



The eigenvalue equation reads 



A 4- AC a 
muo 2 = A + 2C + AD - {3B + 28E) — — — - 4D cos 2 (±) 

-D + lb-C/ 2 

±2(-C + 10E ^ )cos(^). (4.8) 
v B+16E' K 2 J v ' 



Results are summarized in Table |T|. Note that for C = 10 E g^^ E the two branches coincide. 



D. States with period N > 3 

For the N = 3 solution exact phonon frequencies can in principle be found. The elements 
of the dynamical matrix are given in terms of £ and k, defined in ( |3.5| ) and ( |3.6| ): 

Di,i = D 3 ,s = A + 2(C + D) + [3Bk 2 + 2E{Ak 2 - 3k + 1)] £ 2 , 
,02,2 = A + 2{C + D) + [3B + 4£(£; 2 - 3k + 3)] £ 2 , 

£>i, 2 = D 2;i = -C- E(3k 2 -Ak + 3)£ 2 - De~ iq , (4.9) 
D 1>3 = —D — (C + 2Ek 2 i 2 )e~ iq . 

The eigenvalues are then found as the solution (Cardano's formula) of a cubic equation. 
For N = A the dynamical matrix has elements in terms of p, defined in ( |3.9| ): 

D n>n = A + 2(C + D) + (3B + 16£)p 2 , 
£>i,2 = D 3A = -C- 2Ep 2 , 

D lt3 = D 2A = -D(l + e-^), (4.10) 
D 14 = (-C - 10Ep 2 )e- i9 , 
D 2)3 = -C- lOEp 2 . 
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The resulting secular equation is a quartic one and exact solutions for the eigenvalues can 
be found using Ferrari's formula. 

For solutions with larger periods we have to rely on numerical calculations. As an exam- 
ple we again consider the ground state for A = 2.24999, B = 1, (7=1, D = — 1, E = 1. See 
also the end of Section III. The calculated phonon dispersion curves in the commensurate 
approximation A = || are given in Figure |2|. 



V. CALCULATION OF PHASE DIAGRAMS 

In this section we present some phase diagrams, calculated partly analytically, partly 
numerically. The traditional method to find the ground state numerically is to solve the 
equations for equilibrium (|3.1|). However, these equations also hold for metastable states, 
maxima and saddlepoints and it may happen that one finds a metastable state instead of the 
true ground state. This problem is not present for the so-called effective potential method 
(EPM), introduced by Griffiths and Chou.il This method in principle always gives the 
ground state. Originally it was used to study Frenkel-Kontorova and similar one-dimensional 
models with only nearest-neighbor interaction. As an interesting application of this method, 
we mention a study of the ground state of the chiral XY model in a field.il Below we give 
a brief outline of the method. 

Consider a one-dimensional system with only nearest- neighbor interaction in its ground 
state. If one atom is displaced from its equilibrium position (we assume that x n denotes 
the displacement), the surrounding atoms will change their positions in order to minimize 
the total energy. This deformation will in general cost some energy. A function, called the 
'effective potential', describes the net energy cost as a function of the positions of the atoms. 
This effective potential achieves its minimum on points of the ground state" and rigorous 
mathematical statements can be madelll Numerical procedures to find solutions are based 
on discretization of the range x n of the atomic positions. The x n can now only adopt a finite 
number of values. For models with interactions up to next-nearest neighbors, as in the 
case of the extended DIFFOUR model, the EPM can be adapted, which will be discussed in 
Appendix B. The proofs of the existence of solutions for models with next-nearest-neighbor 
interactions, both in the continuous and discretized version, are rather long and will be given 
in a separate paper.E3 

Using the EPM and equation (3.1) we calculated various phase diagrams. First we varied 
both A and E with the other parameters fixed: B — 1, C — 1 and D = — 1. The resulting 
phase diagram is given in Figure |3|. From the analysis in Appendix A we know that the phase 
boundary for the paraelectric state for \C\ < 4 is given by A = jC 2 — 2C + 4. For C = 1 and 
D = — 1 we find A = 2\. At this boundary we have a transition to an incommensurate state 
with wave vector q = arccos(^) = arccos(^) ~ 4.76679213. This is the state we discussed 
at the end of Section III. We can clearly see the effect of the E-term: for E = further 
decreasing of A leads to a transition to a commensurate state with period 4. For E < 
this transition can be followed by transitions to iV = 3 or iV = 2 commensurate states. 
For E sufficiently positive, the wavelength of the ground state increases for decreasing A. 
Between the commensurate states incommensurate ones can be found. By increasing E, 
the region between paraelectric and ferroelectric phases shrinks. A positive E term favors 
long- wavelength solutions, with the ferroelectric state being the extreme limit (q — 0). 
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Figure ^ gives the phase diagram found for B = 1, D = — 1, E = and varying both 
A and C. This is the phase diagram for the original DIFFOUR model.i We have seen that 
the phase boundary for the paraelectric phase for \C\ < 4 is a parabola symmetric around 
C = 0. For \C\ > 4 this boundary is given by A + 2C - 2 = -2 + 2\C\, two straight lines. 
The parabola and the lines meet at |CJ_= 4, and have equal derivative at this point. Note 
the symmetry C <-> — C, which impliest3 that the modulation wave vectors for the system 
with +C and — C are related by 

1 . . 

q c + q- C = - (5.1) 

in units of 2n. At (C = 4, A — 0) the paraelectric phase, the ferroelectric phase and the 
incommensurate phase become equal. The lines separating the paraelectric phase from the 
incommensurate phase and the incommensurate phase from the ferroelectric phase have 
equal derivative at this point. In Landau theory (see Section VII) such a point would be 
called a Lifshitz point. From the symmetry C <-> — C it is obvious that there is also a Lifshitz 
point at (C = —4, A = 16). At this point the paraelectric phase, the antiferroelectric phase 
and the incommensurate phase become equal. 

Figure [| gives the phase diagram for B — 1, D — — 1, E — 1 in terms of A and C. The 
symmetry C <-> —C is no longer present. However, the phase boundary of the paraelectric 
phase is independent of E. Also the wavelength of the phase emanating from this boundary 
is the same. In particular the positions of the Lifshitz points and the derivates at these 
points do not change. Note the boundary of the antiferroelectric phase: starting from the 
Lifshitz point and going down in the phase diagram, it initially bends to the right and then 
returns to lower values of C. Figures £| and [5] have been obtained by solving equation ( |3.1|) 
and comparing the energies of the solutions. 

We investigated the nearby surroundings of the Lifshitz point at (C = 4, A = 0) to look 
how the transition line from the ferroelectric phase to the incommensurate phase changes by 
increasing E. See figure ^| for the results. One notices a tendency towards a smaller wedge 
W (the vertical distance between the paraelectric-incommensurate phase boundary and the 
incommensurate-ferroelectric phase boundary) by increasing E, as was to be expected. 

VI. TEMPERATURE DEPENDENT BEHAVIOR 

As the model under consideration is one- dimensional with short-range interactions, there 
is no phase transition possible at T ^ 0. If we however consider weakly interacting linear 
chains in a three dimensional system, this system can be described by (|2.1| ) as well if we 
interpret the variables x n as averages over planes perpendicular to a fixed direction (the 
c-axis). Phase transitions become possible due to inter-chain couplings. 

To study the temperature dependence of the parameters we take the thermal average of 
the conditions for equilibrium ( [3.1|) , resulting in 

A(x n ) + B(xl) + C (2(x n ) - (x n -i) - (x„+i» + D (2(x n ) - (x n „ 2 ) ~ {x n+2 )) 

+E [A(xD - 3<4z n _i) - 3{x 2 n x n+1 ) + 2{x n x 2 n _ 1 ) + 2{x n x 2 n+1 ) - (x 3 ^) - (x 3 n+1 )] = (6.1) 

We have to distinguish between ground states with {x n } = and ground states with {x n } ^ 
0, where the {x n } are solutions of Q3.1]) . 
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In the former case we assume that the thermal fluctuations of the displacement x n do not 
depend on the lattice site, (x 2 ) — (x n ) 2 ps (x 2 ^) — (x m ) 2 , and if we furthermore approximate 
the correlations by (x 2 x m ) ~ (x 2 )(x m ), the following holds: 

E [— 3(x n x n _i) — 3(x n x n+ i) + 2(x n x n _ 1 ) + 2(x n x n+1 ) — (x n _ 1 ) — (x n+1 )] 
+ (B + AE)(x 3 n ) 

« -E [-3(x n > 2 ((x n _i> + (x n+1 » + 2(x n ) ((x n _x) 2 + (x n+1 > 2 ) - ((x„_i) 3 + (x n+ i) 3 )] 
+4£ ((x 2 > - (x n > 2 ) (2(x n ) - (x n _i) - (x n+1 » 

+(B + AE)(x n ) 3 + B ((xl) - (x n ) 2 ) (x n ). (6.2) 

Inserting the last expression in (|6.1|) we see that the conditions for equilibrium for the 
thermal average of the displacement, (x n ), have the same form as those for the displacements 
x n themselves; the only difference being the replacement of the parameters A and C by 
temperature dependent ones: 

A -> A + BT, 

C + AET, (6.3) 

where T = (x 2 ) — (x n ) 2 is a measure of the thermal fluctuations. So a change in temperature 
will renormalize both parameters A and C (unlike in the DIFFOUR model with E — 0). 

For all other ground state solutions (x„ ^ 0) we calculate the thermal averages around 
x n , where the {x n } satisfy 



; ,-. , J x » x m e e ux w ux m = ( x p\( x i \ fK 4) 

//e-/ 3 (^- i «) 2 /2e-/3(^-^) 2 /2 da;nda;m 

where (3 = 1/T. Three different integrals have to be calculated, yielding 

( X n) = %n + 3x n T 

(x 2 n )=x 2 n + T (6.5) 

Substitution in ( |6.1| ) and comparison with ( j3.1| ) then leads to 

A -> A + (35 + 4£)T, 

C -► C + GET. (6.6) 

Note that the parameter E now also enters in the temperature dependence of A. This linear 
behavior in T, with a kink at the temperature where the transition from the paraelectric 
phase to the incommensurate or the ferroelectric phase takes place, is corroborated by Monte 
Carlo calculations.0il Some of the results!! are shown in Figure [j]. The results ( |6.3j ) and 
( |6.6D are in sharp contrast with the assumptions made in standard Landau theory, to be 
discussed in Section VII, that there is only one temperature dependent parameter, and that 
its behavior above and below the transition temperature is the same. 

It is now straightforward to calculate the temperature dependent ground states and 
stability limits by making substitutions (|6.3|) and (|6.6|). We especially would like to focus on 
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the phonon branches in the paraelectric and ferroelectric phases. Of experimental interest 
is the ratio between the slopes of the soft phonon mode in the ferroelectric and paraelectric 
phase, the so-called i?-parameter: 

d^7dT| ferro 
dc 2 /dT| para 

Self-consistent renormalized phonon theory!! gives the result R = —2. In experiments^ 
however very often R ^ —2 is found. Taking into account the temperature dependence 
given in ( |6.3|) and Q6.6Q , we find using (|4.5|) and Q4.6|) 

= -2(3£ + 4£)-32£ 2 /i?sin 2 (g/2) 

5+ 16£sin 2 (g/2) ' 1 ' ' 

which, for E — 0, gives (at the center of the Brillouin zone) R = —6 instead of R = —2, and 
for E ytz o can take on any value as long as B + 16.E > is satisfied. Molecular dynamics 
simulations on a 3-dimensional 4 lattice by Padlewski et a/.0 show that R = —2 holds 
only for systems with long-range couplings being in the displacive limit. However, Sollich 
et a/.0 showed that this double limit of displaciveness and long range interaction is not 
necessary if the system is displacive enough: they found R = —2 for a system with only 
nearest-neighbor interactions, thereby questioning Padlewski's claim@ of having studied a 
system in the displacive limit. 



VII. COMPARISON WITH THE CONTINUUM THEORY 

The continuum limit of the extended DIFFOUR model leads to a well-known expansion. 
By replacing the differences in the general order parameter x n by differentials in the principal 
order parameter for ferroelectrics, the polarization P(z), 

{x n %n— l) * I ~\ I 



dz J 



(x n -i + x n+1 -2x n ) 2 -> ( -^-2 J , 
and rearranging terms we arrive at the free energy density 

with a = A, (3 = B, k = Cjt AD, X = —D and rj = 2E. The above free energy density was 
used by Ishibashi and ShibaP^ to study phase transitions in NaNC>2 and SC(NH2)2 (thiourea), 
proper ferroelectrics in which the polarization component of interest transforms according to 
a one-dimensional irreducible representation. The r^-term is allowed by symmetry because it 
is the product of the two invariants P 2 and (d 2 P/dz 2 ) 2 . Alternatively, as both sodium nitrite 
and thiourea admit an interaction of P with another mode u (strain for example), Dvofaklil 
showed that the //-term accounts in an effective manner for this interaction, thereby reducing 
g(P,u)^f(P). 
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By taking the Fourier transform of the above free energy density we find 



M^ 2 wHi + ?' 2 ) p « 4 (73) 

This justifies the choice of the P-term in the extended version of the DIFFOUR model, 
discussed in Section II. A term oc (dP/dz) 4 has been included in the free energy expansion 
by Jacobs et a/.,@ but by taking the Fourier transform one finds oc q 4 P 4 which is of higher 
order than the rj-term used here. 

In a seminal paper Hornreich et a/El discussed a multicritical point of a new type, which 
they called a Lifshitz point. In the spherical model limit they were able to calculate critical 
exponents and the shape of the phase boundaries of 2 nd order and 1 st order transitions in 
the vicinity of the Lifshitz point .e^I Let us return to the above free energy density to give a 
definition^ of the Lifshitz point. At an ordinary paraelectric to ferroelectric phase transition 
the coefficient a changes sign. If we have an additional incommensurate phase we need the k 
and A-terms, and at the Lifshitz point k = 0. Higher-order terms in the expansion are needed 
for stabilization. Converting a = 0, k = to variables in the extended DIFFOUR model we 
find A = 0,C = 4 (for D = —1). This is exactly the position of the Lifshitz point found 
in Section V. There is another analogy between Landau theory and the DIFFOUR model: 
MichelsonS showed that for systems with uniaxial polarization the phase transition lines 
separating the paraelectric phase from the incommensurate phase and the incommensurate 
phase from the ferroelectric phase are tangent at the Lifshitz point. This feature is also 
present in figures | and [5]. 

Let us now discuss some properties of the solutions found in Landau theory. Ground 
states minimize the total free energy 

1 r d 



F=- I f(z)dz, (7.4) 



o 



and can be found by solving the Euler-Lagrange equation 

+ aP + (3P 3 = 0. (7.5) 



, d 4 P d 2 P 

A——; K- 



dz 4 dz 2 



fdp\ 2 9 d 2 p 
p — +P 2 



\dz J dz 



2 



GolovkoE3 was able to obtain exact solutions for some special values of the parameters in 
a slightly more general free energy density (he added a term ^P 6 to the expansion ( |7. 2| ) ) . 
However, his method is not general and we will not discuss it further. Instead we follow 
a different approach: numerically solvingil equation ( |7.5| ) shows that the solutions contain 
practically only one harmonic, the amplitude of higher harmonics is at most 3.5% of the 
former. 

As usual in Landau theory only the coefficient a is temperature dependent: a = a Q (T — 
T c ). It is found that between the high-temperature paraelectric solution 

p( z ) =0, F = 0, (7.6) 

and the low-temperature ferroelectric solution 
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/3' " 4/T v "' 

an incommensurate solution exists. Justbelow the paraelectric-incommensurate transition 
at aii = a (Ti — T c ) — it has the formea 

P W =*««<„), F=-^^.. (7.8) 
The amplitude po an d wave vector g are given by 



2 



4(ao — a) 
3/3+2 Wo 2 



2 • 



^o(l + ^Po] (7.9) 



The r^-term makes the incommensurate phase less stable when i] is positive, implying that 
the transition temperature from the incommensurate state to the ferroelectric state increases 
as rj increases. See also the discussion by Toledano.i In the discrete model a positive .E-term 
is responsible for this effect. 



VIII. CONCLUSIONS AND OUTLOOK 

In this paper we have calculated various properties of an extension of the DIFFOUR 
model. For this purpose a next-nearest-neighbor generalization of the Effective Potential 
Method was developed. The shape of the paraelectric phase boundary was proven rigorously, 
elaborating on a former proof which only included nearest-neighbor interactions. We found 
that the phase diagram changes considerably due to the extra E-term, but the transition 
at the paraelectric phase boundary does not depend on E. Positive E favors longer-period 
solutions. 

By taking thermal fluctuations in two different regimes into account the parameters A and 
C can be considered as effectively temperature dependent. For C this holds only for nonzero 
E, which explains the relevance of this extra term. This has strong consequences for two 
experimentally easy accessible quantities: the temperature dependence of the modulation 
wave vector, and the ratio between the slopes of the soft phonon mode in the ferroelectric 
and paraelectric phases (i?-parameter). 

Although lattice and continuum models have some features in common, the differences 
are more striking. A lattice model would be a more natural choice than the phenomenological 
Landau treatment of incommensurate phases. Discrete models do not need ad hoc lock- 
in terms to explain different commensurate and incommensurate phases. Complex phase 
diagrams can in principle be obtained using a simple Hamiltonian which takes into account 
the discreteness of a lattice. 

The Sn2P2(Si_ x Sea;)6 crystal family seems to be an excellent system for our future re- 
search: it is uniaxial, has an exceptional Lifshitz point in the composition-temperature phase 
diagram, shows cross-over effects from order-disorder to a displacive type of phase transition, 
and displays an interesting modulation wave vector behavior. All these phenomena can in 
principle be explained by the extended DIFFOUR model. 
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APPENDIX A: EXACT RESULTS FOR PARAELECTRIC AND 
FERROELECTRIC PHASES 

In this Appendix explicitly calculated phase boundaries are given for the extended DIF- 
FOUR model. We first consider E = and then discuss the effect of E ^ 0. Let us start 
with writing V in the form 

V = S { i x ™ + \ x ™ + cx « x «-i + dx n%n-2 \ , (Al) 



with d = ±1. The remaining parameters a,c,d are the tilde parameters defined in ( |2.2| 
after normalization of B and D. Let us now try to write this as 



V = ^2 { ~ <l x n-i ~ rx n _ 2 ) 2 + -a£ 1 



(A2) 



2 i „2\ a 



Comparison of the two expressions yields 

p(l + q 2 + r / y 

p{-2q + 2qr) = c, (A3) 
— 2pr = d. 

From this one can see that if it is possible to write the potential in this form and a is positive, 
then p is positive. Eliminating q and r from the above equations yields the following fourth 
order polynomial equation (assuming non-zero a, c and d): 

16p 4 + (16d - 8a)p 3 + (8d 2 + Ac 2 - 8ad)p 2 + (4rf 3 - 2ad 2 )p + d A = 0. (A4) 



First consider the d — +1 case. Equation (|A4|) then has the following complex solutions 



V 



V 



U - 2 + a + y/(2 + a) 2 - 4c 2 

±V2^/-4 + a 2 - 2c 2 + (-2 + a) y^T a) 2 - 4c 2 ^ , 

V-2 + a- v/(2 + a) 2 -4c 2 

iv^^-4 + a 2 - 2c 2 - (-2 + a) y^T a) 2 - 4c 2 ^ . 



(A5) 



(A6) 
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The first requirement for having a real solution is that (2 + a) 2 — 4c 2 > 0, i.e. a > —2 + 2\c\. 
Consider the first two solutions (|A5[) . First look at |c| > 4. Then for a > —2 + 2\c\ we find 
(—2 + a)i/(2 + a) 2 — 4c 2 > 0. And for the first term in the root we find — 4 + a 2 — 2c 2 > 
2c 2 — 8|c| > 0. So for \c\ > 4 the only requirement for having a real (positive, a > 0) solution 
is a > -2 + 2|c|. For |c| < 4 we have -4 + a 2 - 2c 2 + (-2 + a)y/(2 + a) 2 - 4c 2 = for 
a = 2 + jc 2 . It can be seen that the argument of the root is positive for a > 2 + ^c 2 . So, 
for |c| < 4 the requirement for having a real (positive, a > 0) solution is: a > 2 + |c 2 (then 
a > — 2 + 2|c| automatically holds too). 

The requirements for |c| > 4 and |c| < 4 form a continuous line in the a — c-parameter 
space. Above this line V can be written as ( |A^ ) with p positive. Therefore V > 0. The 
lower bound is reached by the trivial solution which always exists, so the paraelectric phase 
is the ground state above this line. In Sec. IV it is shown that the above line corresponds 
exactly with the stability lines of the trivial solution, showing that the modulated phases 
arise from the destabilization of the trivial solution due to the condensation of a soft phonon 
mode. 

In case d — — 1 the solutions of the fourth order polynomial equation are 
V = + a + y/(-2 + a) 2 - 4c 2 

±V2\J-A + a 2 - 2c 2 + (2 + a) y/(-2 + a) 2 - 4c 2 ^) , (A7) 



V 



X - {2 + a - y/(-2 + a) 2 - 4c 2 

±V2^/-4 + a 2 - 2c 2 - (2 + a) V / H2~+ a) 2 - 4c 2 ^ . 



(A8) 



Look at the two solutions ( |A7| ). The first condition is a > 2 + 2\c\. (2 + a) is positive if 
this requirement is fulfilled. Further —4 + a 2 — 2c 2 > 2c 2 + 8|c| > 0. So, here we have only 
one requirement for all c, namely a > 2 + 2\c\. Above this line the paraelectric phase is the 
ground state. 

With the same sort of reasoning we can also try to prove that for certain parameter 
values the ground state is ferroelectric Here we work with the following form of V (with 
D = ±1): 

V = H\ ^ x2 n + \ X n + j( X n ~ X n-lf + y (x n ~ X n _ 2 ) 2 1 . (A9) 



We try to write this as 



V = ^ + \ X n + P ( X n - QXn-1 ~ -RX„_ 2 ) 2 | 



(AlO) 



Comparison of the two expressions yields 
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P(1 + Q 2 + R 2 ) = C + D, 

P(-2Q + 2QR) = -C, (All) 
-2PR = -D. 

If there exists a solution and C + D is positive, then P is positive. Rewriting the above 
equations yields the following fourth order polynomial equation (assuming nonzero C and 
D): 

16P 4 + (-32D - 16C)P 3 + (24P 2 + 16CP> + AC 2 )P 2 + (-8P 3 - ACD 2 )P + D A = 0. 

(A12) 

For the case D — — 1 the complex solutions are 

P = i (-2 + C ± + C) , (A13) 

both having multiplicity 2. For < C < A the solution is not real. For C > 4 the solution is 
real. In order to have a positive solution we must have C + D > 0, so C > 1. So, for C > 4 
the potential can be written as ( |A10| ) with P positive. So 

( A14 ) 

The ferroelectric phase, which exists if A < 0, reaches this lower bound. So for B = 1, D = 
— 1, A < the ground state is ferroelectric for C > 4. In terms of the tilde parameters: for 
B — 1, D = 1, A < —2 — 2C, the ferroelectric phase is the ground state for C < —4. 
For D = +1 the complex solutions are 

p = - (2 + C ± v / C\ / 4~Tc) . (A15) 

For —4 < C < the solution is not real. For C > the solution is real. In that case P is 
positive. For P = 1,P = 1,A<0 the ferroelectric phase is the ground state for C > 0. In 
other words: for P = 1,P = — 1, A < 2 — 2C the ferroelectric phase is the ground state for 
C < 0. In terms of the tilde parameters analogous statements about the anti-ferroelectric 
phase can easily be made (C <-> — C). 

In case P 7^ the following holds: the parts of the phase diagram where the paraelectric 
phase is the ground state in the DIFFOUR model (P = 0) also belong to the paraelectric 
phase for this extended model for all allowed values of P. In this extended model there are 
no other parts of the phase diagram where the trivial solution is the ground state, because 
the stability conditions for this solution are the same as in the DIFFOUR model (see Section 
IV). 

For the ferroelectric phase the statements are less rigorous: if the ferroelectric phase is 
the ground state in the DIFFOUR model (P = 0), then it is also the ground state in the 
extended model for P > 0. Again we can prove that V > ^2 n + jx 4 }. However, for 
positive P the part of the phase diagram where the ferroelectric phase is the ground state 
becomes bigger. 
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APPENDIX B: EFFECTIVE POTENTIAL METHOD FOR NEXT-NEAREST 

NEIGHBORS 



In this Appendix, which is based on the account given by Griffiths,0 we discuss how the 
EPM can be generalized to be applicable to systems in which there is next-nearest neighbor 
interaction. Consider a classical one-dimensional chain of atoms. The total potential energy 
of the system is given by 

oo 

H= J2 {V(x n ) + W(x n+1 ,x n ) + D(x n+2 ,x n )} . (Bl) 

n=— oo 

So, interactions up to next-nearest-neighbors are included. The effective potential method 
that is used to find ground states for systems with interaction with first neighbors, can 
also be used for systems where second neighbor interaction is included. Instead of a scalar 
variable at site n. one now has to deal with a vector consisting of the values for x for two 
adjacent atoms.El Writing x n = (x 2n ,x 2n -i) the above potential energy is of the form 

i/ = ^iv(x n ,x n _ 1 ), (B2) 

n 

with 



fT(x n ,x n _i) = V(x 2n ) + V(x 2n -i) + W(x 2 n,X 2 n-l) + W(x 2n -i,x 2n - 2 ) 

+D(x 2n , x 2n ^ 2 ) + D{x 2n ^i, x 2n s). (B3) 

The fact that here a vector at site n is considered does not change the EPM and the proofs 
given for this method-SiH^l The solution can be found by solving 

7] + R(x n ) = min {^(x^i) + K(x n , x^)} . (B4) 

Here 77 is 2 times the ground state energy per particle. The vector consists of two compo- 
nents with respect to which one has to minimize. Because of this minimization over two 
components, which has to be performed frequently, the numerical procedures based on dis- 
cretization of the range of possible x n -values will take a very long time. The method we 
used is slightly different. Consider the n th couple of adjacent atoms. Couple (n — 1) does 
not consist of two other atoms as is the case above. Instead, the right atom of couple (n — 1) 
is the same as the left atom of couple n. In this case only a minimization over one atomic 
degree of freedom (the 'position') is left. This leads to more reasonable computation times. 
The fact that couple (n — 1) is not independent of couple n requires an adaptation of the 
proof of the existence of a solution both in the continuous case and in the discretized case 
used for numerical procedures. We only give an outline for the method, proofs of the ex- 
istence of solutions and generalization to systems with interactions up to s th neighbors are 
to be given in a separate paper.il The method for deriving the equations and the numerical 
procedures^ for solving the equations remain essentially the same. Numerical calculations 
suggest that the error in rj has a cubic dependence on the grid size rather than a quadratic 
dependence for the Frenkel-Kontorova model.il 



18 



Let us give the following explanation^] for the method: Imagine that a system described 
by (|B1|) is in its ground state. If we now change the positions of two adjacent atoms, the 
surrounding atoms will in general also change their positions in order to minimize the total 
energy. This net energy change caused by the deformation of one couple will be called 
the effective two-particle-potential. This will describe the energy cost as a function of the 
positions of two adjacent atoms. At site n, the effective two-particle-potential R(x n+ i,x n ), 
due to the presence of the atoms i < n, can be formally written as 



R(x n+1 ,x n ) =min\ V] Wfa) + W(x h + D(x h x^ 2 ) - v] 

U<n+1 



(B5) 



where the minimum is taken over all atomic positions Xi with i < n and rj is the (unknown) 
ground state energy per particle. By rewriting this equation, one obtains 

R(x n+ i, x n ) = min min <^ S~] [V(xi) + W(x h x^i) + D(x i} x^ 2 ) - v] 

x n -i i<n-l I * — • 
v i<n 

+V{x n+1 ) + W{x n+1 , x n ) + D( (B6) 



which gives 

7] + R(x n+1 , X n ) = V(x n+1 ) + W(x n+ l, X n ) + mm [R(x n , Xn-l) + D(x n +l,X n -l)} • (B7) 



•En. — l 



This is the minimization eigenvalue equation for R. The same procedure can be followed for 
the effect of the atoms i > n + 1. The effective two-particle-potential due to these atoms is 
called S(x n +i,x n ), which gives 

r] + S(x n+ i, x n ) = V(x n ) + W(x n+1 , x n ) + min [S(x n+2 , x n+ i) + D(x n +2, x n )] . (B8) 

The total effective two-particle-potential F(x n +i,x n ) of a couple of adjacent atoms in a 
double infinite chain, is given by 

) ~ V(x n+ i) - V(x n ) - W(x n+ i, x n ), (B9) 

where the last three terms are subtracted on the right side to avoid double counting. 

The equations (|B7|) and ( |B8|) can also be obtained in another wayHH Let R,N(x n +i,x n ) 



be the minimal energy of a chain of iV atoms with the constraint that the atoms iV and 
N — 1 are at fixed positions x n+ \ and x n respectively, while the other atoms are free to 
rearrange themselves in an optimal way so as to minimize the total energy. This leads to 

R 2 {x n+ 1, X n ) = V{x n+ l) + V{x n ) + W(x n+1 , X n ), (BIO) 
) = V(x n +l) + V(x n ) + W(x n+ l, X n ) 

+ min \V{x n „i) + W{x n , x n -i) + D(x n +i, x n ~i)] 
= V(x n+ i) + W(x n+1 ,x n ) + min R 2 {x n , x n _i) + D( 

), (Bll) 

%n — 1 L 

) = V(x n+ i) + W(x n+1 , x n ) + min R N {x n , x n ^) + D(x n+ i, x n -i) ■ (B12) 
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Now assume that for iV — > oo, R N (x n+ i, x n ) approaches some function R(x n+ i,x n ) plus a 
constant proportional to iV — 2: 

R N (x n+1 , x n ) -> R(x n+1 ,x n ) + {N - 2)rj. (B13) 



In that case equation ( |B7|) follows. However, it is not clear that ( B13 ) will always be satisfied. 
But by imposing a special boundary condition, namely 

R 2 (x n+1 , x n ) = R(x n+1 , x n ), (B14) 



( P13|) will be satisfied exactly.^ The previous boundary condition is the same as saying 



that the left-most couple experiences the effective two-particle-potential instead of the true 
two-particle-potential. The minimum energy of this system as a function of the positions 
of the two right-most atoms is given by R + Nrj. Assuming that R is a bounded function, 
the energy per particle of such a system will tend to 77 as iV — > 00. 77 is thus the average 
energy per particle in any ground state, since the extra boundary condition only changes 
the total energy by a term of order 1. So R is the effective two-particle-potential for the 
right-most couple of a semi-infinite chain. The same is true for S for the left-most couple 
of a semi-infinite chain extending to the right. F is the total effective two-particle-potential 
for a couple in a double-infinite chain. R,S and F can of course only be defined up to an 
additive constant. 

In the above derivations the problems arising from the summation of an infinite number 
of terms in ( |B5| ) are neglected. In fact, one considers local deformations of length M, with 
the limit M — > 00. This will be explained below. Define the effective two-particle-potential 
due to the local deformation of length M as 

R iM) (x n+1 ,x n ) = min <^ V] [K{x h x^ u x^ 2 ) - v] 

n+l—M<i<n I ' 

K n+3-M<i<n+l 

+ [K[X n+3 -M, x n+2-M, W„ + i_Af) 

+K(x n+2 -M,Un+l-M,U n - M ) ~ 2T]] J , ( B1 5) 

where Ui refers to the ground state value for atom i, and where we have introduced 

K(x n+1 ,x n , x n -i) = V(x n+1 ) + W(x n+1 ,x n ) + D(x n +i,x n -i). (B16) 
The right hand site of equation ( |B15| ) can be rewritten as 

R {M \x n+1 ,x n ) = min [i? (M_1) (x„, x n -i) + K(x n+1 ,x n ,x n ^ 1 ) - 77] . (B17) 

It is reasonable to assume that in the limit M — > 00 : R^ M \x n+ i, x n ) — > 

R (M-1) ( 

(because x n+2 -M -> u h+ 2-m)- Writing R(x n+1 ,x n ) = lim A f_oo R {M) {x n+1 , x n ) the mini- 
mization eigenvalue equation results. In the second version of obtaining the equations it is 
clear that it is in fact the limit of local deformations, however with the boundary condition 
that the left most couple of atoms experiences the effective two-particle-potential. Here, 
one should take the length of the chain going to infinity in order to let rj go to the ground 
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state energy per particle. The above explanation also holds for S. It is best to picture the 
situation as a local deformation of the ground state. 

Now, the nonlinear minimization eigenvalue equations for R and S are rewritten.!!!! 
The eigenvalue equation for R now becomes 

77 + R(x n+1 ,x n ) = min [R(x n , x n _i) + K(x n+1 ,x n , x n ~i)] ■ (B18) 

Let the function L be defined by 

L(x n+1 , x n ) = S(x n+1 , x n ) - V(x n+ i) - V(x n ) - W(x n+ i, x n ). (B19) 
The minimization eigenvalue equation for S can now be rewritten as 

77 + L(x n+1 ,x n ) = min [L(x n+2 , x n+1 ) + K(x n+2 , x n+1 , x n )} . (B20) 

Xn + 2 

In terms of R and L one has 

Ffen+l) x n) = R\ x n+li x n) ~t~ L[X n +i, X n ). (B21) 

In fact there may be multiple solutions of the eigenvalue equation, not only differing 
by a trivial constant. The existence of different solutions is related to the existence of 
different degenerate ground states. The general solution is given by 

R(x n+ i,x n ) = min [R a (x n+1 ,x n ) + K a ] . (B22) 

a 

The R a correspond to the pure phases and the K a are arbitrary constants. 

For each solution of the minimization eigenvalue equation for R ( |B18j ), a r map can 



be defined,S0 where r(x n+ i,x n ) = {(x n ,x n _i)} with x n _i one of the values for which the 
minimum on the right hand side of ( |B18| ) is achieved. An i?-orbit is defined as 

Vn : {x n , x n -i) G T{x n +x, x n ) t] + R{x n +i, x n ) = R{x n , x n -\) + K(x n +\, x n , x n _i). 

(B23) 



Similarly, for the minimization eigenvalue equation for L (|B2(Jp, a o map can be defined, 



where a(x n+ i,x n ) = {(x n+ 2, x n+ i)} with x n+ 2 one of the values for which the minimum on 
the right hand side of flB20f) is achieved. An L-orbit is defined as 



Vn : (x n+2 , x n+ i) e cr(x n+1 , x n ) =>• i] + L(x n+1 , x ri ) = L(x n+2 , x n+1 ) + K(x n+2 , x n+1 , x n ). 

(B24) 

A ground state is both an i?-orbit and an L-orbit. Therefore, it can be proven that for a 
ground state 

F(x n+ i, x n ) = R(x n+1 ,x n ) + L(x n+ i,x n ) = F(x n , x n _i) (B25) 

So, F is constant on the positions of two adjacent atoms in a ground state, which is logical 
since it is the effective two-particle-potential. 
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Numerical procedures are based on a discretized version of the system. In that case for 
each ground state there is a solution for the eigenvalue equation for which there is a path 
from each point to the ground state in the corresponding r graph. El So, the situation is as 
follows. There is a local deformation of length M (with M — > oo), in a chain coinciding 
with a particular ground state for ±00. This ground state corresponds to a certain solution 
of the eigenvalue equation. The deformation is such that atoms n and n + 1 have values x n 
and x n+ \. By applying the corresponding r map one can obtain the positions of the atoms 
left from n. For the discretized system one will finally reach the ground state in this way 
(In the continuous case it is supposed to converge to the ground state). The same can be 
done for the atoms right from n + 1 by applying the a map. From this picture it is clear 
that the ground state is both an R- orbit and an L-orbit. In fact the r map and the a 
map may be multi-valued. So, the atomic positions of the atoms (say) left from n do not 
have to be unique. The deformation can have parts consisting of minimizing cycles (cycles 
of minimal energy) different from the ground state configuration at ±00. In the r graph 
one can go directly to the minimizing cycle corresponding to the ground state configuration 
at —00, or one can first stay for some time in another minimizing cycle if this exists. If 
there are several solutions for R and S (with several corresponding r and a maps), there are 
several possibilities to construct F. It will often be logical to take the ground states, toward 
which the chain converges at —00 and +00, the same. The numerical algorith we used is an 
adapted version of the one discussed by Floria and Griffiths^ 

Here an example will be given to show that it is important that in fact limits of finite 
deformations are considered. Suppose that the ground state is ferroelectric, with two de- 
generate ground states: x n = x = ±/ where 1^0. If the deformation is just infinite as 
suggested in ( [Bop the value of R(l, I) should be the same as the value for R(—l, —I). How- 
ever, something else is seen. Two solutions can be found corresponding to the two ground 
states. In the solution corresponding to the solution x n = I, R(—l, —I) has a higher value 
than R(l, I). The difference is the defect energy, the energy cost for going from the +1 phase 
to the — I phase. (The defect energy (and the defect configuration) can also be calculated 
using the r map.) From this it can be seen that the deformation is in fact embedded in the 
ground state Ui = +1 at —00 (In the limit M — ► 00: x n+2 -M U n+2~M where u n+2 -M = +1, 
or for the second version: the left most couple of atoms in the finite chain experiences the 
effective two-particle-potential corresponding to the ground state Ui = +1). 

It can also be expected that F has local minima at the positions of two adjacent atoms 
in metastable states. However, since only two atoms are at a fixed position, while the other 
atoms are free to rearrange themselves in an optimal way, this may not be the case. When 
changing the positions of two adjacent atoms by an infinitesimal amount, the changes of the 
other atomic positions in the metastable state does not have to be infinitesimal. Therefore, 
it is not necessarily true that there is a local minimum in F for positions of two adjacent 
atoms in a metastable state. When there are no other atomic positions in the metastable 
state (a period 1 solution), F does have a local minimum. When the lowest metastable state 
has positions of two adjacent atoms which are not seen in a ground state (which will often 
be the case), there will be a local minimum in F for these two positions. In that case the 
energy cannot be lowered by changing the other atoms by any amount, since the only states 
which have lower energy are ground states and these cannot be reached since the positions 
of the two adjacent atoms in consideration are not in a ground state (and the changes of 
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them should be infinitesimal) . By following the development of the shape of F one may also 
investigate the kind of phase transitions that are involved. For example, a discontinuous 
change in the set of points where F achieves its global minimum, indicates a first order 
transition.^ 
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FIGURES 




FIG. 1. Ground state configuration for A = 2.24999, B = 1, C = 1, L> = -1 and £ = 1 in 
the commensurate approximation 62/13. x n is the displacement for particle n. For A, B,C, D, E 
given above the system is just below the phase boundary between the paraelectric phase and 
the incommensurate phase. Although the displacements are small, 0(1O -3 ), the onset of the 
incommensurate phase is evident. Points are calculations, lines are drawn to guide the eye. 
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FIG. 2. Phonon dispersion curves for A = 2.24999, B = 1, C = 1, D = -1 and E = 1 in the 
commensurate approximation 62/13. This corresponds with the solution depicted in Figure |l[ q is 
given in reduced units. There is one branch with oj — > for q — > 0: the phason branch. Just above 
this lies the amplitudon branch. 
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FIG. 3. The phase diagram for the extended DIFFOUR model with B = 1, C = 1, D = -1. 
Shown are the paraelectric phase (P), ferroelectric phase (F), antiferroelectric phase (2), and com- 
mensurate phases with period 3,4,5,6. Incommensurate phases are labeled I, higher order commen- 
surate phases C. There is no stable ground state for E < —1/16. 
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FIG. 4. The phase diagram for the extended DIFFOUR model with B = 1, D = -1 and E = 0. 
This is the original DIFFOUR model. Note the symmetry C «-* — C. Same labeling as in figure ||. 
L denotes the Lifshitz points. 
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FIG. 5. The phase diagram for the extended DIFFOUR model with B = 1, D = — 1 and E = 1. 
Note the asymmetric character, although the boundary of the paraelectric phase is the same as in 
figure Same labeling as in figure |I[ 
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FIG. 7. Temperature dependence A(T) calculated on a 3D version of (2.1) with harmonic 



nearest-neighbor coupling and A(0) = —5, B = 5, D = E = 0. Points are calculations,! 
fits with a linear function. The phase transition takes place at A = 0. 



lines are 



32 



TABLES 

TABLE I. Stability limits for the paraelectric phase ( |4.5| ). Note that B + 16E > must be 

satisfied. 

Parameter range g-value of instability Conditions for having a stable state 

D < 0: 

C < AD q c = TT A> -AC 

AD<C<-AD cos( 9c ) = =g A+ (*c+i6D) 2 >Q 

C > -AD q c = A > 

D > 0: 

C < g c = vr ,4 > -4(7 

C >0 g c = A > 



TABLE II. Stability limits for the ferroelectric phase 
must be satisfied. 



. Note that B + IQE > and A < 



Parameter range 



g-value of instability 



Conditions for having a stable state 



D < 0: 



C< 



2AE 
B 



+ AD 



2AE + 4jD < C < 2AE 
B 



c > 



AD 



AD cos(q c ) - ^jj -r 



C _|_ AE 



A < 2C - ^ 

0A , (C+AD—2AE/B) 2 

A < 



> 



I? > 0: 

r< / 2AE 
U < -g- 

^ ^ 2A£ 



g c = o 



A < 2C - ^ 
A < 



TABLE III. Stability limits for the antiferroelectric phase 
A + AC < must be satisfied. 



II). Both B + 16E > and 



Parameter range 



q- value of instability 



Conditions for having a stable state 



D < 0: 
C ~ WE-fcffig < AD 



q c = 



A < -AC 



AD<C - WE^±^ < cos(f 



< C - WE^±^ < -AD cos( 



IZJ-^iWM) A>-2C-AD + (3B + 28E)^ 

A+4C 



Ml 1 ""' y 1/J 



4£> \^ w T B+16E ^ 

A > -2C -AD + (3B + 28£)ig|| 



C 



1QT7i^4+4C_ ^ 4 J") 



L> > 0: 



A+4C 



< 



B+16E 
3B+8E 
B+16E 



A < - 



A+4C 



-AC 

3B+8E 
B+16E 
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